Introduction

https://github.com/edwardlavender/patter

Introduction

Team

Institutions

Workshop

• The goal today is to introduce patter
• We will focus on how to use the package
• There will be:
  • Interactive presentation & questions
  • (Time allowing) Optional practical activity (using example datasets or your own)

Questionaire

• Experience with animal-tracking?
• Experience with R and Julia?
• Experience with statistical modelling, including SSMs?
• Experience with Bayesian statistics?
• Experience with particle algorithms?

Context

• patter uses a dedicated high-performance backend written in Julia

https://github.com/edwardlavender/Patter.jl

Evolution

• patter evolved from flapper (Lavender et al., 2023)

Evolution

• flapper is a more general-purpose package for passive acoustic telemetry that supports:
  • data acquisition, simulation, assembly and processing
  • spatial operations and distance metrics (e.g., fast C++ LCP algorithms)
  • data exploration (detection statistics, movement metrics)
  • ‘flapper’ algorithms

Evolution

• patter focuses specifically on model inference for SSMs using particle algorithms

https://github.com/edwardlavender/patter

Methodology

1. We formulate a Bayesian state-space model for the joint probability distribution \(f(\mathbf{s}_{1:T} \mid \mathbf{y}_{1:T})\) of a tagged individual’s possible trajectories from the start to the end of the time series

2. The joint probability distribution is expressed in terms of a prior (the movement process) and a likelihood (the observation process)

\[f(\mathbf{s}_{1:T} \mid \mathbf{y}_{1:T}) \propto f(\mathbf{s}_{1:T}) f(\mathbf{y}_{1:T} \mid \mathbf{s}_{1:T})\]

Methodology

3. Model inference is performed using particle filtering and smoothing

4. The particle filter approximates the partial marginal distribution \(f(\mathbf{s}_t \mid \mathbf{y}_{1:t})\) using a set of weighted samples called ‘particles’

5. The particle smoother approximates the full marginal distribution \(f(\mathbf{s}_t \mid \mathbf{y}_{1:T})\)

Methodology

Inference using particle algorithms.

Methodology

Particle filter intuition.

Methodology

• Particle smoothers reweight particles from the filter to approximate \(f(\mathbf{s}_{t} \mid \mathbf{y}_{1:T})\)
• The two-filter smoother consists of a re-weighting of particles from a forward filter run and a backward filter run based on the probability density of movements between pairs of locations.

Methodology

The two filter smoother.

Operating systems

• patter v.1.0.0 supports Windows and MacOS
• patter@dev supports Linux and will be rolled out soon

Installation

• To install patter, follow the instructions on https://github.com/edwardlavender/patter
  • Install R and Julia
  • On Windows, install RTools
  • On MacOS, install Xcode (and configure .Makevars as required)
  • Install patter (and let it handle Julia dependencies)

Connect to Julia

• The syntax for julia_connect() is:

julia_connect(..., JULIA_PROJ, .threads, ...)

• In patter@dev, this is:

julia_connect(JULIA_HOME, JULIA_PROJ, JULIA_NUM_THREADS, ...)

• If you’ve set environmental variables (recommended), just use:

julia_connect()

Overview

## Overview

Overview

• We will work though:
  • Boilerplate setup
  • An example workflow with real data
  • An example simulation-based workflow

Overview

patter workflow

Boilerplate setup

• Connect to Julia:

julia_connect()

• Set the seed in R and Julia via set_seed():

set_seed(2024)

Study system

Observations include detections at receivers:

head(dat_acoustics)

   individual_id           timestamp receiver_id
           <int>              <POSc>       <int>
1:            25 2016-03-17 01:50:00          26
2:            25 2016-03-17 01:52:00          26
3:            25 2016-03-17 01:54:00          26
4:            25 2016-03-17 01:58:00          26
5:            25 2016-03-17 02:00:00          26
6:            25 2016-03-17 02:04:00          26

Study system

We also have archival (depth) observations, collected at a resolution of two-minutes:

head(dat_archival)

   individual_id           timestamp  depth
           <num>              <POSc>  <num>
1:            25 2016-03-16 00:00:00 159.15
2:            25 2016-03-16 00:02:00 161.01
3:            25 2016-03-16 00:04:00 162.86
4:            25 2016-03-16 00:06:00 164.48
5:            25 2016-03-16 00:08:00 158.46
6:            25 2016-03-16 00:10:00 155.45

Study system

• Our goal in the first part of this presentation is to reconstruct movements and patterns of space use for an example individual given these data.

# Define detections
det <-
  dat_acoustics |>
  filter(individual_id == 25L) |>
  mutate(individual_id = NULL) |>
  as.data.table()

# Define archival (depth) observations
arc <-
  dat_archival |>
  filter(individual_id == 25L) |>
  mutate(individual_id = NULL) |>
  as.data.table()

Study system

plot(arc$timestamp, arc$depth * -1, 
     col = "royalblue", type = "l",
     ylim = c(-max(arc$depth), 0), 
     xlab = "Time stamp", ylab = "Depth (m)")
points(det$timestamp, rep(0, nrow(det)))

Study system

• We need to define a map of our study system in R and Julia

Study system

Study system

• Cell entries should be floats
• Code inhospitable habitats (e.g., land) as NA_real_
• A planar (UTM) coordinate reference system is required

Study system

• Define the study timeline:
  • This sets the period, and resolution, at which we represent animal movement
  • We may have observations regularly or irregularly along this timeline

Study system

• Use assemble_timeline() to set a timeline based on the input datasets:

timeline <- assemble_timeline(.datasets = list(det, arc), 
                              .step = "2 mins", 
                              .trim = TRUE)
str(timeline)

 POSIXct[1:24225], format: "2016-03-17 01:50:00" "2016-03-17 01:52:00" "2016-03-17 01:54:00" ...

Study system

• Set the timeline using seq.POSIXt():

timeline <- seq(as.POSIXct("2016-03-17 01:50:00", tz = "UTC"),
                as.POSIXct("2016-03-18 01:48:00", tz = "UTC"), 
                by = "2 mins")
str(timeline)

 POSIXct[1:720], format: "2016-03-17 01:50:00" "2016-03-17 01:52:00" "2016-03-17 01:54:00" ...

State

• State is an Abstract Type in Patter.jl
• This is an object that holds the parameters describing the ‘state’ of an individual
• Typically, state means ‘location’ (in two or three dimensions), but additional parameters may be included

# The Julia code for the StateXY struct
struct StateXY <: State
    map_value::Float64
    x::Float64
    y::Float64
end 

State

• In Patter.jl, the following sub-types are built-in:

  • StateXY(map_value, x, y): Used for two dimensional (x, y) states ;
  • StateXYZD(map_value, x, y, z, angle): Used for four-dimensional (x, y, z, direction) states;

State

• In patter, you just need to specify the state sub-type as a character string:
  • "StateXY" maps onto StateXY in Julia
  • "StateXYZD" maps onto StateXYDZ in Julia

state <- "StateXY"

• Each state-type maps onto a corresponding movement model

Movement model

• ModelMove is an Abstract Type in Patter.jl
• This is an object that holds the dimensions of the movement model, e.g.,
  • The map
  • A distribution of step lengths
  • A distribution of turning angles

# The Julia code for the ModelMoveXY struct
struct ModelMoveXY{T, U, V} <: ModelMove
    map::T
    dbn_length::U
    dbn_angle::V
end 

Movement model

• The movement model can be instantiated via an R move_*() wrapper
  • StateXY maps onto ModelMoveXY, which is wrapped by move_xy()
  • StateXYZD maps onto ModelMoveXYZD, which is wrapped by move_xyzd()
• The fields of the movement model in the move_*() wrapper must be specified using Julia code

Movement model

• Define a two-dimensional random walk (x, y):

# Define maximum moveable speed between two time steps
mobility   <- 750.0

# Define 2D RW (XY)
model_move <- move_xy(dbn_length = glue("truncated(Gamma(1.0, 250.0), upper = {mobility})"),
                      dbn_angle = "Uniform(-pi, pi)")

Movement model

• move_*() wrappers return a string of Julia code:

model_move

ModelMoveXY(env, truncated(Gamma(1.0, 250.0), upper = 750), Uniform(-pi, pi));

Movement model

• Visualise the components of our movement model:

expr <- expression({
  pp <- par(mfrow = c(1, 2))
  curve(dtrunc(x, spec = "gamma", a = 0, b = mobility, shape = 1.0, scale = 250.0), 
        from = 0, to = mobility, n = 1e3L,
        xlab = "Step length (m)", ylab = "Density")
  curve(dunif(x, -pi, pi),
        from = -pi - 0.1, to = pi + 0.1, n = 1e3L,
        xlab = "Heading (rad)", ylab = "Density")
  par(pp)
})

Movement model

• Visualise the components of our movement model:

Movement model

• Visualise realisations of the movement model:

Movement model

• More sophisticated movement models are possible, including:
  • Random walks in two or three dimensions
  • Correlated random walks
  • Behavioural-switching models
• See the help for ?State to get started

Movement model

• This is an example of a more sophisticated movement model for flapper skate:

Example behavioural-switching model.

Observation model(s)

• We have defined the movement model
• This encodes the process of animal movement
• We now need to define the observation models
• These encode the probability of the observations, given the latent (unobserved) location of the animal

Observation model(s)

• As a reminder, we have collected detections and depth time series:

head(det, 2)

             timestamp receiver_id
                <POSc>       <int>
1: 2016-03-17 01:50:00          26
2: 2016-03-17 01:52:00          26

head(arc, 2)

             timestamp  depth
                <POSc>  <num>
1: 2016-03-16 00:00:00 159.15
2: 2016-03-16 00:02:00 161.01

Observation model(s)

• Each dataset should comprise a time series of observations for one individual
• Datasets should be pre-processed as required: e.g.,
  • Check receiver locations are valid on the map
  • Identify false detections
• patter accepts any kind of observational dataset

Observation model(s)

• Each dataset should be provided as a data.table with the following columns:
  • sensor_id, timestamp, obs;
  • Additional columns that define observation-model parameters;

Observation model(s)

• In passive acoustic telemetry datasets, we record detections but the observations comprise detections and non-detections
• Use the assemble_acoustics() helper to assemble the acoustic (detection, non-detection) time series:

acc <- assemble_acoustics(.timeline  = timeline,
                          .acoustics = det,
                          .moorings  = dat_moorings, 
                          .services  = NULL)

Observation model(s)

• If you have passive acoustic telemetry observations, you should also use acoustic containers
• The syntax is:

# ModelObsAcousticContainer
containers <- assemble_acoustics_containers(.acoustics, .mobility, ...)

Observation model(s)

• We have to assemble other datasets with the required columns (sensor_id, timestamp, obs, ...) manually

head(arc, 2)

             timestamp  depth
                <POSc>  <num>
1: 2016-03-16 00:00:00 159.15
2: 2016-03-16 00:02:00 161.01

arc <- 
  arc |> 
  # Define 'dummy' sensor_id
  mutate(sensor_id = 1L) |> 
  # Observations must be in the 'obs' column
  rename(obs = "depth") |> 
  # Select columns
  select("sensor_id", "timestamp", "obs") |>
  # Add observation model parameters 
  mutate(depth_sigma = 50, depth_deep_eps = 20) |>
  as.data.table()

Observation model(s)

• It is a good idea to pass other datasets through the assemble_archival() function for processing/checks though

arc <- assemble_archival(.timeline = timeline, .archival = arc)
head(arc)

   sensor_id           timestamp   obs depth_sigma depth_deep_eps
       <int>              <POSc> <num>       <num>          <num>
1:         1 2016-03-17 01:50:00 73.78          50             20
2:         1 2016-03-17 01:52:00 73.32          50             20
3:         1 2016-03-17 01:54:00 73.32          50             20
4:         1 2016-03-17 01:56:00 73.32          50             20
5:         1 2016-03-17 01:58:00 73.55          50             20
6:         1 2016-03-17 02:00:00 68.70          50             20

Observation model(s)

• For each dataset, we need to define the corresponding ModelObs sub-type
• A ModelObs sub-type is a structure that holds the parameters of an observation model
• ModelObs sub-types are instantiated in Julia using the data.table(s) assembled above
• We use ModelObs instances to simulate observations or to evaluate the log-probability of observations, given latent location

Observation model(s)

• The following ModelObs structures are built-in to Patter.jl:
  • ModelObsAcousticLogisTrunc
  • ModelObsDepthNormalTrunc
  • ModelObsDepthUniform

Observation model(s)

• In a truncated logistic detection probability model:
  • receiver_alpha and receiver_beta affect the rate of decline in detection probability with distance from a receiver
  • receiver_gamma is the maximum detection range

Observation model(s)

Example truncated logistic detection probability models

Observation model(s)

• You can explore different parameters values using this code:

expr <- expression({
  # Define parameters 
  receiver_alpha <- 4
  receiver_beta  <- -0.01
  receiver_gamma <- 1000
  # Visualise model 
  curve(ifelse(x <= receiver_gamma, 
              1 / (1 + exp(-(receiver_alpha + receiver_beta * x))), 
              0), 
       from = 0, to = receiver_gamma, 
       xlab = "Distance (m)", ylab = "Detection probability")
})

Observation model(s)

eval(expr)

Observation model(s)

• ModelObsDepthNormalTrunc is a ModelObs structure for a depth observation and a truncated normal model
• This is a suitable depth-observation model for benthic/demersal species
• The individual must be located in an envelope around the bathymetric depth, defined by a normal distribution centred at this location, truncated at depth_deep_eps m below the seabed and 0 m, with standard deviation depth_sigma

Observation model(s)

Example truncated logistic detection probability models

Observation model(s)

• ModelObsDepthUniform is ModelObs structure for a depth observation and a uniform depth model
• This model assumes that an individual must be located in an envelope around the bathymetric depth, defined by two error terms (depth_shallow_eps and depth_deep_eps)
• This model can be tailored for benthic or pelagic species

Observation model(s)

• For the particle filter, multiple sets of observations should be provided as a named list
• List names define the model sub-types

yobs <- list(ModelObsAcousticLogisTrunc = acc, 
             ModelObsDepthNormalTrunc = arc)

Observation model(s)

• As for the movement model, the parameters of the observation model(s) are assumed known

Particle filter

• We have now defined:
  1. Study area map
  2. Study timeline
  3. State
  4. Movement model
  5. Observations/observation model parameters/observation model sub-types
• We are now in a position to run the particle filter!

Two-filter smoother

• Run the two-filter smoother using outputs from pf_filter() in the Julia environment:

smo <- pf_smoother_two_filter(.n_particle = 100L, .n_sim = 100L)

Two-filter smoother

• pf_smoother_two_filter() returns the familiar pf_particles-class object:

class(smo)

[1] "list"         "pf_particles"

summary(smo)

            Length Class      Mode   
xinit       0      -none-     NULL   
states      6      data.table list   
diagnostics 4      data.table list   
convergence 1      -none-     logical

Two-filter smoother

• We can analyse the outputs in the same way as those from the filter:

head(smo$states)

   path_id timestep           timestamp map_value        x       y
     <int>    <int>              <POSc>     <num>    <num>   <num>
1:       1        1 2016-03-17 01:50:00  73.27666 709144.9 6253203
2:       1        2 2016-03-17 01:52:00  91.54701 709011.6 6253335
3:       1        3 2016-03-17 01:54:00  73.27666 709165.8 6253162
4:       1        4 2016-03-17 01:56:00  56.53801 709377.6 6253224
5:       1        5 2016-03-17 01:58:00  73.44244 709305.7 6253476
6:       1        6 2016-03-17 02:00:00  91.54701 709088.4 6253267

Two-filter smoother

• diagnostics includes ESS:

head(smo$diagnostics)

   timestep           timestamp      ess maxlp
      <int>              <POSc>    <num> <num>
1:        1 2016-03-17 01:50:00      NaN   NaN
2:        2 2016-03-17 01:52:00 76.60220   NaN
3:        3 2016-03-17 01:54:00 78.19307   NaN
4:        4 2016-03-17 01:56:00 50.58620   NaN
5:        5 2016-03-17 01:58:00 81.76317   NaN
6:        6 2016-03-17 02:00:00 81.01046   NaN

Two-filter smoother

• convergence is necessarily TRUE:

smo$convergence

[1] TRUE

Mapping

• We can compute probability-of-use (POU) using particle samples from the filter (for comparison) or the smoother (recommended):

map_pou(.map = map, .coord = smo$states)$ud

Mapping

Mapping

• In practice, POU is sensitive to grid resolution and kernel smoothing may be desirable

map_dens(
  .map,
  .owin = as.owin.SpatRaster(.map),
  .coord = NULL,
  .discretise = FALSE,
  .shortcut = list(),
  ...,
  .fterra = FALSE,
  .plot = TRUE,
  .use_tryCatch = TRUE,
  .verbose = getOption("patter.verbose")
)

Mapping

Residency

• We can easily compute residency in a particular area (at any arbitrary time-scale of our choosing) as the proportion of particles inside/outside an area:

# Define an example polygon (e.g., MPA)
spatpoly <- 
  cbind(708433.9, 6251765) |> 
  terra::vect(crs = terra::crs(map)) |> 
  terra::buffer(width = 2000)

# Visualise study area and 'MPA'
terra::plot(map)
points(smo$states$x, smo$states$y, pch = ".", col = "red")
terra::lines(spatpoly, col = "black", lwd = 2)

Residency

Residency

# Lookup whether or not particles are in the polygon
spatcoord <- terra::vect(cbind(smo$states$x, smo$states$y))
smo$states[, in_poly := 
               terra::relate(spatcoord, spatpoly, relation = "intersects")]
               
# Calculate the expected time spent in the 'MPA' (%)
table(smo$states$in_poly) / nrow(smo$states) * 100

   FALSE     TRUE 
71.89583 28.10417 

Boilerplate setup (1/7)

# Load & attach package 
library(patter)

# Load datasets
head(det, 2)
head(arc, 2)

# Connect to Julia 
julia_connect()
set_seed()

# Set map 
map <- dat_gebco()
set_map()

# Study timeline 
timeline <- assemble_timeline(.datasets = list(det, arc))

Set state (2/7)

state <- "StateXY"

Set movement model (3/7)

mobility   <- 750.0
model_move <- move_xy(dbn_length = glue("truncated(Gamma(1.0, 250.0), upper = {mobility})"),
                      dbn_angle = "Uniform(-pi, pi)")

Set observation model(s) (4/7)

# Acoustic observations 
acc <- assemble_acoustics(.timeline  = timeline,
                          .acoustics = det,
                          .moorings  = dat_moorings, 
                          .services  = NULL)

# Archival observations 
arc <- assemble_archival(.timeline = timeline, 
                         .archival = 
                           arc |> 
                             rename(obs = "depth") |> 
                             select("sensor_id", "timestamp", "obs") |>
                             mutate(depth_sigma = 50, depth_deep_eps = 20) |> 
                           as.data.table())

# All observations
yobs <- list(ModelObsAcousticLogisTrunc = acc, 
             ModelObsDepthNormalTrunc = arc)

Run particle filter (5/7)

# List filter arguments
args <- list(.map         = map,
             .timeline    = timeline,
             .state       = state,
             .xinit_pars  = list(mobility = mobility),
             .yobs        = yobs,
             .model_obs   = names(yobs),
             .model_move  = model_move,
             .n_particle  = 1e5L)

# Forward run
fwd <- do.call(pf_filter, args, quote = TRUE)

# Backward run
args$.direction <- "backward"
bwd <- do.call(pf_filter, args, quote = TRUE)

Run particle smoother (6/7)

smo <- pf_smoother_two_filter()

Analyse outputs (7/7)

map_dens(.map = map, .coord = smo, sigma = bw.diggle)

Simulation-based workflow

• patter also supports simulation-based workflows:
  1. sim_array() simulates arrays
  2. sim_path_walk() simulates movement paths
  3. sim_observations() simulates observations
• You can model simulated observations

Simulate arrays

• Use sim_array() to simulate randomly or regularly arranged receivers
• Locations are sampled from .map
• Default observation model parameters (receiver_alpha, receiver_beta, receiver_gamma) can be assigned

moorings <- sim_array(.map = map, 
                      .timeline = timeline,
                      .arrangement = "random",
                      .n_receiver = 100L, 
                      .receiver_alpha = 5, 
                      .receiver_beta = -0.01, 
                      .receiver_gamma = 1200, 
                      .n_array = 1L, 
                      .plot = TRUE)

Simulate arrays

Simulate arrays

head(moorings)

   array_id receiver_id      receiver_start        receiver_end receiver_x
      <int>       <int>              <POSc>              <POSc>      <num>
1:        1           1 2016-03-17 01:50:00 2016-03-18 01:48:00   708542.1
2:        1           2 2016-03-17 01:50:00 2016-03-18 01:48:00   705742.1
3:        1           3 2016-03-17 01:50:00 2016-03-18 01:48:00   707842.1
4:        1           4 2016-03-17 01:50:00 2016-03-18 01:48:00   710642.1
5:        1           5 2016-03-17 01:50:00 2016-03-18 01:48:00   706542.1
6:        1           6 2016-03-17 01:50:00 2016-03-18 01:48:00   706642.1
   receiver_y receiver_alpha receiver_beta receiver_gamma
        <num>          <num>         <num>          <num>
1:    6265607              5         -0.01           1200
2:    6266007              5         -0.01           1200
3:    6253007              5         -0.01           1200
4:    6254707              5         -0.01           1200
5:    6255407              5         -0.01           1200
6:    6255707              5         -0.01           1200

Simulate paths

• We have already seen how to simulate paths

# Set state and movement model 
state      <- "StateXY"
mobility   <- 500.0
model_move <- 
  move_xy(dbn_length = "truncated(Normal(50, 250.0), lower = 0, upper = {mobility})", 
                      dbn_angle = "VonMises(0, 0.25)")

# Simulate path using movement model 
pp <- par(mfrow = c(1, 2))
path_coord <- sim_path_walk(.map = map, 
                            .timeline = timeline, 
                            .state = state, 
                            .n_path = 1L, 
                            .model_move = model_move)

# Fit UD using cross validation
path_ud <- map_dens(.map = map, .coord = path_coord, .discretise = TRUE)$ud
par(pp)

Simulate paths

Simulate paths

head(path_coord)

   path_id timestep           timestamp map_value        x       y
     <int>    <int>              <POSc>     <num>    <num>   <num>
1:       1        1 2016-03-17 01:50:00  9.677233 701642.1 6265607
2:       1        2 2016-03-17 01:52:00 14.596510 701931.3 6265519
3:       1        3 2016-03-17 01:54:00 20.246822 702173.1 6265576
4:       1        4 2016-03-17 01:56:00 20.431116 702256.8 6265467
5:       1        5 2016-03-17 01:58:00 20.431116 702260.3 6265465
6:       1        6 2016-03-17 02:00:00 20.246822 702174.9 6265569

Simulate observations

• To simulate observations, we need to choose our observation models and their parameters
• For each observation, we need to provide parameters as a data.table that defines, for each sensor_id (e.g., receiver), the observation-model parameters associated with that sensor

Simulate observations

• Let’s set things up to simulate acoustic observations:

# ModelObsAcousticLogisTrunc 
# > We have multiple acoustic sensors
# > This data.table includes the relevant parameters
pars_ModelObsAcousticLogisTrunc <-
    moorings |>
    select(sensor_id = "receiver_id",
           "receiver_x", "receiver_y",
           "receiver_alpha", "receiver_beta", "receiver_gamma") |>
    as.data.table()

Simulate observations

• And depth observations:

# ModelObsDepthNormalTrunc
# > We have one depth sensory 
# > This data.table includes the relevant parameters 
pars_ModelObsDepthNormalTrunc <- data.frame(sensor_id = 1L,
                                            depth_sigma = 20,
                                            depth_deep_eps = 20)

Simulate observations

• Simulate observations along .timeline via sim_observations():

obs <- sim_observations(
  .timeline = timeline, 
  .model_obs = c("ModelObsAcousticLogisTrunc",
                 "ModelObsDepthNormalTrunc"),
  .model_obs_pars = list(pars_ModelObsAcousticLogisTrunc, 
                         pars_ModelObsDepthNormalTrunc))

Simulate observations

• sim_observations() returns a list, with one element for every .model_obs
• Each element is a list, with one sub-element for each simulated path
• Each sub-element is a data.table that contains the observations

Simulate observations

str(obs)

List of 2
 $ ModelObsAcousticLogisTrunc:List of 1
  ..$ :Classes 'data.table' and 'data.frame':   72000 obs. of  8 variables:
  .. ..$ timestamp     : POSIXct[1:72000], format: "2016-03-17 01:50:00" "2016-03-17 01:50:00" ...
  .. ..$ obs           : int [1:72000] 0 0 0 0 0 0 0 0 0 0 ...
  .. ..$ sensor_id     : int [1:72000] 1 2 3 4 5 6 7 8 9 10 ...
  .. ..$ receiver_x    : num [1:72000] 708542 705742 707842 710642 706542 ...
  .. ..$ receiver_y    : num [1:72000] 6265607 6266007 6253007 6254707 6255407 ...
  .. ..$ receiver_alpha: num [1:72000] 5 5 5 5 5 5 5 5 5 5 ...
  .. ..$ receiver_beta : num [1:72000] -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 -0.01 ...
  .. ..$ receiver_gamma: num [1:72000] 1200 1200 1200 1200 1200 1200 1200 1200 1200 1200 ...
  .. ..- attr(*, ".internal.selfref")=<externalptr> 
 $ ModelObsDepthNormalTrunc  :List of 1
  ..$ :Classes 'data.table' and 'data.frame':   720 obs. of  5 variables:
  .. ..$ timestamp     : POSIXct[1:720], format: "2016-03-17 01:50:00" "2016-03-17 01:52:00" ...
  .. ..$ obs           : num [1:720] 17.8 17.8 13.1 32.3 16 ...
  .. ..$ sensor_id     : int [1:720] 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ depth_sigma   : num [1:720] 20 20 20 20 20 20 20 20 20 20 ...
  .. ..$ depth_deep_eps: num [1:720] 20 20 20 20 20 20 20 20 20 20 ...
  .. ..- attr(*, ".internal.selfref")=<externalptr> 

Simulate observations

acc <- obs$ModelObsAcousticLogisTrunc[[1]]
head(acc)

             timestamp   obs sensor_id receiver_x receiver_y receiver_alpha
                <POSc> <int>     <int>      <num>      <num>          <num>
1: 2016-03-17 01:50:00     0         1   708542.1    6265607              5
2: 2016-03-17 01:50:00     0         2   705742.1    6266007              5
3: 2016-03-17 01:50:00     0         3   707842.1    6253007              5
4: 2016-03-17 01:50:00     0         4   710642.1    6254707              5
5: 2016-03-17 01:50:00     0         5   706542.1    6255407              5
6: 2016-03-17 01:50:00     0         6   706642.1    6255707              5
   receiver_beta receiver_gamma
           <num>          <num>
1:         -0.01           1200
2:         -0.01           1200
3:         -0.01           1200
4:         -0.01           1200
5:         -0.01           1200
6:         -0.01           1200

Simulate observations

arc <- obs$ModelObsDepthNormalTrunc[[1]]
head(arc)

             timestamp      obs sensor_id depth_sigma depth_deep_eps
                <POSc>    <num>     <int>       <num>          <num>
1: 2016-03-17 01:50:00 17.83227         1          20             20
2: 2016-03-17 01:52:00 17.81938         1          20             20
3: 2016-03-17 01:54:00 13.13591         1          20             20
4: 2016-03-17 01:56:00 32.34808         1          20             20
5: 2016-03-17 01:58:00 15.98056         1          20             20
6: 2016-03-17 02:00:00 16.76636         1          20             20

Simulate observations

# Plot depths
  plot(arc$timestamp, arc$obs * -1, 
       col = "royalblue", type = "l",
       ylim = c(-max(arc$obs), 0), 
       xlab = "Time stamp", ylab = "Depth (m)")
# Add detections 
det <- acc[obs == 1L, ]
points(det$timestamp, rep(0, nrow(det)), col = "red")

Simulate observations

Mapping

pp <- par(mfrow = c(1, 3))

# Plot the UD for the simulated path 
terra::plot(path_ud)
map_hr_home(path_ud, .add = TRUE)

# Plot the estimated UD 
patter_ud <- map_dens(.map = map, 
                   .coord = smo$states, 
                   .discretise = TRUE, 
                   .fterra = TRUE,
                   sigma = bw.diggle)$ud
map_hr_home(patter_ud, .add = TRUE)

# Compare to a UD based on the COA algorithm
coa_ud <- map_dens(.map = map, 
                   .coord = coa(.map = map, 
                                .acoustics = det, 
                                .moorings = moorings, 
                                .delta_t = "2 hours", 
                                .plot_weights = FALSE), 
                   .discretise = FALSE, 
                   .fterra = TRUE, 
                   sigma = bw.diggle)$ud

par(pp)

Mapping

Custom State, ModelMove ModelObs

• To incorporate custom structures, associated simulation and density Julia methods are required:
  • Patter.simulate_step() and Patter.logpdf_step()
  • Patter.simulate_obs() and Patter.logpdf_obs()

Custom State and ModelMove

• For custom State and ModelMove structures, we need to:
  1. Register the structures in Julia
  2. (optional) Write states_init() methods to translate initial coordinates to States
  3. Write simulation and density methods:
     • Patter.simulate_step(state, model, ...) (for simulation)
     • logpdf_step(state_from, state_to, model_move, ...) (for smoothing)
  4. (optional) Write a move_*() wrapper for convenient R usage

Custom ModelObs

• For custom ModelObs structures, we need to:
  1. (optional) Write map_init() method to sample initial States efficiently
  2. Write simulation and density methods:
     • Patter.simulate_obs(state, model, ...)
     • Patter.logpdf_obs(state, model, ...)

Custom workflow

• Define a custom State
• We will track the location of an animal in three dimensions:

struct StateXYZ <: Patter.State
  # Map value
  # > This is required for all states & is the value on the map at (x, y)
  map_value::Float64
  # Coordinates
  x::Float64
  y::Float64
  z::Float64
end

Custom workflow

• (optional) Define a function to convert initial coordinates to the state
• Here, we simply have to add initial depth values
• This is necessary for automated sampling of initial states

 states_init.StateXYZ <- function(.state, .coords) {
  # Define global variables
  z <- map_value <- NULL
  # Add initial z values
  .coords[, z := map_value * runif(.N)]
  .coords
}

Custom workflow

• Define a corresponding movement model

struct ModelMoveXYZ{T, U, V, W} <: Patter.ModelMove
  # The environment (i.e., map)
  # > This defines the regions within which movements are permitted (i.e., in water)
  map::T
  # Distribution for step lengths
  dbn_length::U
  # Distribution for turning angles
  dbn_angle::V
  # Distribution for changes in depth
  dbn_z_delta::W
end

Custom workflow

• Define a function to simulate new states using the movement model
• This is necessary for sim_path_walk() and pf_filter()

function Patter.simulate_step(state::StateXYZ, model::ModelMoveXYZ, t::Int64)
  # Simulate a step length
  length = rand(model.dbn_length)
  # Simulate a turning angle
  angle  = rand(model.dbn_angle)
  # Calculate new x, y, z coordinates
  x      = state.x + length * cos(angle)
  y      = state.y + length * sin(angle)
  z      = state.z + rand(model.dbn_z_delta)
  # Include the map value
  map_value = extract(model.map, x, y)
  # Define the state
  StateXYZ(map_value, x, y, z)
end

Custom workflow

• Provide a logpdf_step() method to compute the unnormalised log-probability of an unrestricted move
• This is necessary for pf_filter_two_filter()

function Patter.logpdf_step(state_from::StateXYZ, state_to::StateXYZ,
                            model_move::ModelMoveXYZ,
                            t::Int64,
                            length::Float64, angle::Float64)
  # Calculate change in depth
  z_delta = state_to.z - state_from.z
  # Evaluate unnormalised log probability
  logpdf(model_move.dbn_length, length) +
    logpdf(model_move.dbn_angle, angle) +
    logpdf(model_move.dbn_z_delta, z_delta)
end

Custom workflow

• (optional) Write an R wrapper to conveniently trial different parameterisations from R:

move_xyz <- function(length = "truncated(Gamma(1.0, 750.0), upper = 750.0)",
                     angle = "Uniform(-pi, pi)",
                     z_delta = "Normal(0, 3.5)") {
  glue::glue('ModelMoveXYZ(env, {length}, {angle}, {z_delta});')
}

Custom workflow

• Now we’ll define a custom observation model for our pelagic species
• We imagine we have depth observations

struct ModelObsDepthNormal <: Patter.ModelObs
  sensor_id::Int64
  depth_sigma::Float64
end

Custom workflow

• (optional) Write a map_init() method for more efficient automated sampling of initial states
• This example routine ensures that we only sample from locations where the bathymetric depth is at least as deep at the observed depth

map_init.ModelObsDepthNormal <- 
  function(.map, .timeline, .direction, .dataset, .model, .pars) {
    t1 <- ifelse(.direction == "forward", .timeline[1], .timeline[length(.timeline)])
    pos <- which(.dataset$timestamp == t1)
    if (length(pos) == 0L) {
      return(.map)
    }
    depth <- .dataset$obs[pos]
    terra::mask(.map, .map >= depth, maskvalues = 0)
}

Custom workflow

• (optional) Write a Patter.simulate_obs() method to simulate new observations
• This is necessary for sim_observations()

function Patter.simulate_obs(state::StateXYZ, model::ModelObsDepthNormal, t::Int64)
  dbn   = truncated(Normal(state.z, model.depth_sigma), 0, state.map_value)
  rand(dbn)
end

Custom workflow

• Write a Patter.logpdf_obs() method to evaluate the log-probability of the observations
• This is required for pf_filter()

function Patter.logpdf_obs(state::State, model::ModelObsDepthNormal, t::Int64, obs::Float64)
  dbn   = truncated(Normal(state.map_value, model.depth_sigma),
                    0.0, state.map_value)
  logpdf(dbn, obs)
end

Custom workflow

• Simulate an example movement path:

 state      <- "StateXYZ"
 mobility   <- 750.0
 model_move <- move_xyz()
 path <- sim_path_walk(.map = map,
                       .timeline = timeline,
                       .state = state,
                       .model_move = model_move)

Custom workflow

• Simulate an example movement path:

Custom workflow

• Simulate observations:

# Run simulation 
 obs <-
   sim_observations(.timeline =
                      timeline,
                    .model_obs =
                      "ModelObsDepthNormal",
                    .model_obs_pars =
                      list(data.table(sensor_id = 1L, depth_sigma = 5)))
 arc  <- obs$ModelObsDepthNormal[[1]]
 yobs <- list(ModelObsDepthNormal = arc)

 # Visualise simulated observations 
 plot(arc$timestamp, arc$obs * -1,
      col = "royalblue", type = "l",
      ylim = c(-max(arc$obs), 0),
      xlab = "Time stamp", ylab = "Depth (m)")

Custom workflow

Custom workflow

• Run the algorithms:

 args <- list(.map         = map,
              .timeline    = timeline,
              .state       = state,
              .xinit_pars  = list(mobility = mobility),
              .yobs        = yobs,
              .model_obs   = names(yobs),
              .model_move  = model_move,
              .n_particle  = 1e5L)

 fwd <- do.call(pf_filter, args, quote = TRUE)
 args$.direction <- "backward"
 bwd <- do.call(pf_filter, args, quote = TRUE)
 smo <- pf_smoother_two_filter()

Custom workflow

• Compare simulated depths with particle depths:

plot(smo$states$timestamp, smo$states$z * -1, 
     pch = ".", col = "dimgrey", 
     ylim = c(-max(smo$states$z), 0),
     xlab = "Time stamp", ylab = "Depth (m)")
lines(arc$timestamp, arc$obs * -1, col = "royalblue")

Custom workflow

• Compare simulated and reconstructed patterns of space use:

pp <- par(mfrow = c(1, 2))
map_dens(.map = map, 
         .coord = path, 
         .discretise = TRUE, 
         sigma = bw.diggle)
map_dens(.map = map, 
         .coord = smo$states, 
         .discretise = TRUE, 
         sigma = bw.diggle)
par(pp)

Custom workflow

Parallelisation

• In practice, we want to analyse the data for many individuals

• You have two options for parallelisation:

  1. (Recommended) Iterate over individuals in sequence and use multi-threading in Julia
  2. Parallelise over time series in R in a socket cluster and use single-threading in Julia

• Forking is not supported with JuliaCall

Parallelisation

• For large-scale analyses, it is a good idea to gauge (and minimise) requirements
• These rack up quickly if you are not careful